home *** CD-ROM | disk | FTP | other *** search
/ SGI Hot Mix 17 / Hot Mix 17.iso / HM17_SGI / research / lib / recon3.pro < prev    next >
Encoding:
Text File  |  1997-07-08  |  17.9 KB  |  434 lines
  1. ; $Id: recon3.pro,v 1.11 1997/01/30 23:48:45 dan Exp $
  2. ;
  3. ; Copyright (c) 1993-1997, Research Systems, Inc.  All rights reserved.
  4. ;       Unauthorized reproduction prohibited.
  5. ;
  6. ;+
  7. ; NAME:
  8. ;    RECON3
  9. ;
  10. ; PURPOSE:
  11. ;    This function can reconstruct a 3-dimensional data array from
  12. ;       two or more images (or projections) of an object.   For example,
  13. ;       if you placed a dark object in front of a white background and
  14. ;       then photographed it three times (each time rotating the object a
  15. ;       known amount) then these three images could be used with RECON3
  16. ;       to approximate a 3-D volumetric representation of the object.
  17. ;       RECON3 also works with translucent projections of an object.
  18. ;       RECON3 returns a 3-D byte array.   RECON3 uses the back-projection
  19. ;       method.   In medical imaging and other applications, a method
  20. ;       known as "Filtered Backprojection" is often desired.   This may
  21. ;       be accomplished here by first filtering the images as desired,
  22. ;       and then using the filtered images for the reconstruction.
  23. ;
  24. ; CATEGORY:
  25. ;    Volume Reconstruction
  26. ;
  27. ; CALLING SEQUENCE:
  28. ;       vol = RECON3(Images, Obj_Rot, Obj_Pos, Focal, Dist, $
  29. ;                    Vol_Pos, Img_Ref, Img_Mag, Vol_Size)
  30. ;
  31. ; INPUTS:
  32. ;       Images:   The images to use to reconstruct the volume.   Execution
  33. ;                 time increases linearly with more images.
  34. ;                 Data Type: 8-bit (byte) array with dimensions [x, y, n]
  35. ;                 where x is the horizontal image dimension, y is the vertical
  36. ;                 image dimension, and n is the number of images.
  37. ;
  38. ;       Obj_Rot:  The the amount the object is rotated to make it appear as
  39. ;                 it does in each image.   The object is first rotated
  40. ;                 about the X axis, then about the Y axis, and finally
  41. ;                 about the Z axis (with the object's reference point at the
  42. ;                 origin.
  43. ;                 Data Type: [3, n] Float array where Obj_Rot[0, *] is the X
  44. ;                 rotation for each image, Obj_Rot[1, *] is the Y rotation, 
  45. ;                 and Obj_Rot[2, *] is the Z rotation.
  46. ;                 
  47. ;       Obj_Pos:  The position of the the object's reference point RELATIVE to
  48. ;                 the camera lens.   The camera lens is located at the
  49. ;                 coordinate origin and points in the negative Z direction
  50. ;                 (the view up vector points in the positive Y direction).
  51. ;                 Obj_Pos should be expressed in this coordinate system.
  52. ;                 The values for Obj_Pos, Focal, Dist, and Vol_Pos should all
  53. ;                 be expressed in the same units (mm, cm, m, in, ft, etc.).
  54. ;                 Data Type: [3, n] Float array where Obj_Pos[0, *] is the X
  55. ;                 position for each image, Obj_Pos[1, *] is the Y position, 
  56. ;                 and Obj_Pos[2, *] is the Z position.   All the values in
  57. ;                 Obj_Pos[2, *] should be less than zero.
  58. ;
  59. ;       Focal:    The focal length of the lens for each image.   Focal may be
  60. ;                 set to zero to indicate a parallel image projection
  61. ;                 (infinite focal length).
  62. ;                 Data Type: Float array with n elements.
  63. ;
  64. ;       Dist:     The distance from the camera lens to the image plane (film)
  65. ;                 for each image.   Dist should be greater than Focal.
  66. ;                 Data Type: Float array with n elements.
  67. ;
  68. ;       Vol_Pos:  The two opposite corners of a cube that surrounds the object.
  69. ;                 Vol_Pos should be expressed in the object's coordinate system
  70. ;                 RELATIVE to the object's reference point.
  71. ;                 Data Type: [3, 2] Float array where Vol_Pos[*, 0] specifies
  72. ;                 one corner and Vol_Pos[*, 1] specifies the opposite corner.
  73. ;
  74. ;       Img_Ref:  The pixel location at which the object's reference point
  75. ;                 appears in each of the images.
  76. ;                 Data Type: [2, n] Int or Float array where Img_Ref[0, *] is
  77. ;                 the X coordinate for each image and Img_Ref[1, *] is the Y
  78. ;                 coordinate.
  79. ;
  80. ;       Img_Mag:  The magnification factor for each image.   This number is
  81. ;                 actually the length (in pixels) that a test object would
  82. ;                 appear in an image if it were N units long and N units
  83. ;                 distant from the camera lens.
  84. ;                 Data Type: [2, n] Int or float array where Img_Mag[0, *] is
  85. ;                 the X dimension (in pixels) of a test object for each image,
  86. ;                 and Img_Mag[1, *] is the Y dimension.   All elements in
  87. ;                 Img_Mag should be greater than or equal to 1.
  88. ;
  89. ;       Vol_Size: The size of the volume to return.   The returned volume will
  90. ;                 be a 3-D byte array with dimensions equal to Vol_Size.
  91. ;                 Execution time (and resolution) increases exponentially with
  92. ;                 larger values for Vol_Size.
  93. ;                 Data Type: Int array with 3 elements where Vol_Size[0]
  94. ;                 specifies the X dimension of the volume, Vol_Size[1] specifies
  95. ;                 the Y dimension, and Vol_Size[2] specifies the Z dimension.
  96. ;
  97. ; KEYWORD PARAMETERS:
  98. ;       CUBIC:    If set, then cubic interpolation is used.   The default is
  99. ;                 to use tri-linear interpolation, which is slightly faster.
  100. ;
  101. ;       MISSING:  The value for cells in the 3-D volume that do not map to
  102. ;                 any of the supplied images.   The Missing value is passed
  103. ;                 to the IDL "INTERPOLATE" function.
  104. ;                 Data Type: Byte.
  105. ;                 Default : 0B
  106. ;
  107. ;       MODE:     If Mode is less than zero then each cell in the 3-D volume
  108. ;                 is the MINIMUM of the corresponding pixels in the images.
  109. ;                 If Mode is greater than zero then each cell in the 3-D volume
  110. ;                 is the MAXIMUM of the corresponding pixels in the images.
  111. ;                 If Mode is equal to zero then each cell in the 3-D volume
  112. ;                 is the AVERAGE of the corresponding pixels in the images.
  113. ;                 Mode should usually be (-1) when the images contain a bright
  114. ;                 object in front of a dark background.   Mode should usually
  115. ;                 be (+1) when the images contain a dark object in front of a
  116. ;                 light background.   AVERAGE mode requires more memory since
  117. ;                 the volume array must temporarily be kept as an INT array
  118. ;                 instead of a BYTE array.
  119. ;                 Data Type: Int
  120. ;                 Default : 0 (average cells)
  121. ;
  122. ; OUTPUTS:
  123. ;    RECON3 returns a 3-D byte array containing the reconstructed object.
  124. ;
  125. ;       If the images contain low (dark) values where the object is and high
  126. ;       (bright) values where the object isn't, then Mode should be set to (+1).
  127. ;       If the above is true then the returned volume will have low values
  128. ;       where the object is, and high values where the object isn't.
  129. ;
  130. ;       If the images contain high (bright) values where the object is and low
  131. ;       (dark) values where the object isn't, then Mode should be set to (-1).
  132. ;       If the above is true then the returned volume will have high values
  133. ;       where the object is, and low values where the object isn't.
  134. ;
  135. ; RESTRICTIONS:
  136. ;       In general, the object must be CONVEX for a good reconstruction to be
  137. ;       possible.   Concave regions are not easily reconstructed.
  138. ;       An empty coffee cup, for example, would be reconstructed as if it
  139. ;       were full.
  140. ;
  141. ;    The images should show strong light/dark contrast between the object
  142. ;       and the background.
  143. ;
  144. ;       The more images the better.   Images from many different angles will
  145. ;       improve the quality of the reconstruction.   It is also important to
  146. ;       supply images that are parallel and perpendicular to any axes of
  147. ;       symmetry.   Using the coffee cup as an example, at least one image
  148. ;       should be looking through the opening in the handle.
  149. ;
  150. ;       Telephoto images are also better for reconstruction purposes than
  151. ;       wide angle images.
  152. ;
  153. ; PROCEDURE:
  154. ;    A 4x4 transformation matrix is created for each image based upon the
  155. ;       parameters Obj_Rot, Obj_Pos, Focal, Dist, and Img_Ref.   Each cell in
  156. ;       the volume is assigned a 3-D coordinate based upon the parameters
  157. ;       Vol_Pos and Vol_Size.   These coordinates are multiplied by the
  158. ;       transformation matricies to produce x,y image coordinates.   Each cell
  159. ;       in the volume is assigned a value that is the AVERAGE, MINIMUM, or
  160. ;       MAXIMUM of the image values at the x,y position (depending on Mode).
  161. ;
  162. ; EXAMPLE:
  163. ; ------------------------------------------------------------------------------
  164. ;       ; Assumptions for this example :
  165. ;       ; The object's major axis is parallel to the Z axis.
  166. ;       ; The object's reference point is at its center.
  167. ;       ; The camera lens is pointed directly at this reference point.
  168. ;       ; The reference point is 5000 mm in front of the camera lens.
  169. ;       ; The focal length of the camera lens is 200 mm.
  170. ;
  171. ;       ; If the camera is focused on the reference point, then the
  172. ;       ; distance from the lens to the camera's image plane must be
  173. ;       ;    dist = (d * f) / (d - f) =
  174. ;       ;    (5000 * 200) / (5000 - 200) = (1000000 / 4800) = 208.333 mm
  175. ;
  176. ;       ; The object is roughly 600 mm wide and 600 mm high.
  177. ;       ; The reference point appears in the exact center of each image.
  178. ;
  179. ;       ;  If the object is 600 mm high and 5000 mm distant from the camera
  180. ;       ;  lens, then the object image height must be
  181. ;       ;     hi = (h * f) / (d - f)  =
  182. ;       ;     (600 * 200) / (5000 - 200) = (120000 / 4800) = 25.0 mm
  183. ;       ;  The object image appears 200 pixels high so the final magnification
  184. ;       ;  factor is
  185. ;       ;     img_mag = (200 / 25) = 8.0
  186. ; 
  187. ;
  188. ;       imgy = 256
  189. ;       frames = 3
  190. ; 
  191. ;       images = Bytarr(imgx, imgy, frames, /Nozero)
  192. ;       obj_rot = Fltarr(3, frames)
  193. ;       obj_pos = Fltarr(3, frames)
  194. ;       focal = Fltarr(frames)
  195. ;       dist = Fltarr(frames)
  196. ;       vol_pos = Fltarr(3, 2)
  197. ;       img_ref = Fltarr(2, frames)
  198. ;       img_mag = Fltarr(2, frames)
  199. ; 
  200. ;       vol_size = [40, 40, 40]
  201. ; 
  202. ;       ; The object is 5000 mm directly in front of the camera.
  203. ;       obj_pos[0, *] =     0.0
  204. ;       obj_pos[1, *] =     0.0
  205. ;       obj_pos[2, *] = -5000.0
  206. ; 
  207. ;       ; The focal length of the lens is constant for all the images.
  208. ;       focal[*] = 200.0
  209. ; 
  210. ;       ; The distance from the lens to the image plane is also constant.
  211. ;       dist[*] = 208.333
  212. ; 
  213. ;       ; The cube surrounding the object is 600 mm X 600 mm.
  214. ;       vol_pos[*, 0] = [-300.0, -300.0, -300.0]
  215. ;       vol_pos[*, 1] = [ 300.0,  300.0,  300.0]
  216. ; 
  217. ;       ; The image reference point appears at the center of all the images.
  218. ;       img_ref[0, *] = imgx / 2
  219. ;       img_ref[1, *] = imgy / 2
  220. ; 
  221. ;       ; The image magnification factor is constant for all images.
  222. ;       ; (The images haven't been cropped or resized).
  223. ;       img_mag[*, *] = 8.0
  224. ; 
  225. ;       ; Only the object rotation changes from one image to the next.
  226. ;       ; Note that the object is rotated about the X axis first, then Y,
  227. ;       ; and then Z.
  228. ; 
  229. ;       ; Create some fake images for this example.
  230. ;       images[30:160, 20:230, 0] = 255
  231. ;       images[110:180, 160:180, 0] = 180
  232. ;       obj_rot[*, 0] = [-90.0, 0.0, 0.0]
  233. ; 
  234. ;       images[70:140, 100:130, 1] = 255
  235. ;       obj_rot[*, 1] = [-70.0, 75.0, 0.0]
  236. ; 
  237. ;       images[10:140, 70:170, 2] = 255
  238. ;       images[80:90, 170:240, 2] = 150
  239. ;       obj_rot[*, 1] = [-130.0, 215.0, 0.0]
  240. ; 
  241. ;       ; Reconstruct the volume.
  242. ;       vol = RECON3(images, obj_rot, obj_pos, focal, dist, vol_pos, img_ref, $
  243. ;                    img_mag, vol_size, Missing=255B, Mode=(-1))
  244. ; ------------------------------------------------------------------------------
  245. ;
  246. ; MODIFICATION HISTORY:
  247. ;     Written by:    Daniel Carr    Thu Feb  4 02:54:29 MST 1993
  248. ;    KDB - 23 Dec., 1993 - Variable dist had a conflict with Userlib
  249. ;                  function DIST and could cause compile errors.
  250. ;                  Renamed variable dist to distance.
  251. ;       Modified by:    Daniel Carr     Mon Nov 21 14:21:57 MST 1994
  252. ;          Improved performance and added CUBIC keyword.
  253. ;       Modified by:    Daniel Carr     Tue Nov 22 12:18:15 MST 1994
  254. ;          Fixed bug which affected small focal length images.
  255. ;          Improved performance again.
  256. ;-
  257.  
  258. FUNCTION RECON3, images, obj_rot, obj_pos, focal, distance, vol_pos, img_ref, $
  259.                  img_mag, vol_size, Missing=missing, Mode=mode, Cubic=cubic
  260.  
  261. ; *** Check inputs.
  262.  
  263. sz_images = Size(images)
  264. IF (sz_images[0] NE 3) THEN $
  265.    Message, 'Image array must have 3 dimensions.'
  266. frames = sz_images[3]
  267. ysize = sz_images[2]
  268. xsize = sz_images[1]
  269.  
  270. sz_obj_rot = Size(obj_rot)
  271. IF (sz_obj_rot[0] NE 2) THEN $
  272.    Message, 'Obj_Rot must be a [3, n] array.'
  273. IF (sz_obj_rot[1] NE 3) THEN $
  274.    Message, 'Obj_Rot must be a [3, n] array.'
  275. IF (sz_obj_rot[2] NE frames) THEN $
  276.    Message, 'Obj_Rot must be a [3, n] array, where n is the number of images.'
  277. obj_rot1 = Float(obj_rot)
  278.  
  279. sz_obj_pos = Size(obj_pos)
  280. IF (sz_obj_pos[0] NE 2) THEN $
  281.    Message, 'Obj_Pos must be a [3, n] array.'
  282. IF (sz_obj_pos[1] NE 3) THEN $
  283.    Message, 'Obj_Pos must be a [3, n] array.'
  284. IF (sz_obj_pos[2] NE frames) THEN $
  285.    Message, 'Obj_Pos must be a [3, n] array, where n is the number of images.'
  286. ind = Where(obj_pos[2, *] GE 0.0)
  287. IF (ind[0] GE 0) THEN $
  288.    Message, 'The object Z position must be < 0 for all images.'
  289. obj_pos1 = Float(obj_pos)
  290.  
  291. IF (N_Elements(focal) NE frames) THEN $
  292.    Message, $
  293.       'Focal must contain the same number of elements as there are images.'
  294. ind = Where(focal[*] LT 0.0)
  295. IF (ind[0] GE 0) THEN $
  296.    Message, 'Focal must be >= 0 for all images.'
  297. focal1 = Float(focal[*])
  298.  
  299. IF (N_Elements(distance) NE frames) THEN $
  300.    Message, 'Len must contain the same number of elements as Focal.'
  301. ind = Where(distance[*] LE focal[*])
  302. IF (ind[0] GE 0) THEN $
  303.    Message, 'Len must be greater than Focal for all images.'
  304. dist1 = Float(distance[*])
  305.  
  306. sz_vol_pos = Size(vol_pos)
  307. IF (sz_vol_pos[0] NE 2) THEN $
  308.    Message, 'Vol_Pos must be a [3, 2] array.'
  309. IF (sz_vol_pos[1] NE 3) THEN $
  310.    Message, 'Vol_Pos must be a [3, 2] array.'
  311. IF (sz_vol_pos[2] NE 2) THEN $
  312.    Message, 'Vol_Pos must be a [3, 2] array.'
  313. vol_pos1 = Float(vol_pos)
  314. IF (vol_pos[0, 0] EQ vol_pos[0, 1]) THEN $
  315.    Message, 'Vol_Pos contains invalid X coordinates.'
  316. IF (vol_pos[1, 0] EQ vol_pos[1, 1]) THEN $
  317.    Message, 'Vol_Pos contains invalid Y coordinates.'
  318. IF (vol_pos[2, 0] EQ vol_pos[2, 1]) THEN $
  319.    Message, 'Vol_Pos contains invalid Z coordinates.'
  320. vpx1 = vol_pos[0, 0] < vol_pos[0, 1]
  321. vpy1 = vol_pos[1, 0] < vol_pos[1, 1]
  322. vpz1 = vol_pos[2, 0] < vol_pos[2, 1]
  323. vpx2 = vol_pos[0, 0] > vol_pos[0, 1]
  324. vpy2 = vol_pos[1, 0] > vol_pos[1, 1]
  325. vpz2 = vol_pos[2, 0] > vol_pos[2, 1]
  326.  
  327. sz_img_ref = Size(img_ref)
  328. IF (sz_img_ref[0] NE 2) THEN $
  329.    Message, 'Img_Ref must be a [2, n] array.'
  330. IF (sz_img_ref[1] NE 2) THEN $
  331.    Message, 'Img_Ref must be a [2, n] array.'
  332. IF (sz_img_ref[2] NE frames) THEN $
  333.    Message, 'Img_Ref must be a [2, n] array, where n is the number of images.'
  334. img_ref1 = Float(img_ref)
  335.  
  336. sz_img_mag = Size(img_mag)
  337. IF (sz_img_mag[0] NE 2) THEN $
  338.    Message, 'Img_Mag must be a [2, n] array.'
  339. IF (sz_img_mag[1] NE 2) THEN $
  340.    Message, 'Img_Mag must be a [2, n] array.'
  341. IF (sz_img_mag[2] NE frames) THEN $
  342.    Message, 'Img_Mag must be a [2, n] array, where n is the number of images.'
  343. ind = Where(img_mag[*] LT 1)
  344. IF (ind[0] GE 0L) THEN $
  345.    Message, 'All elements in Img_Mag must be >= 1.'
  346. img_mag1 = Float(img_mag)
  347.  
  348. IF (N_Elements(vol_size) NE 3) THEN $
  349.    Message, 'Vol_size must contain 3 elements.'
  350. vol_size1 = Long(Abs(vol_size[*]))
  351.  
  352. ; *** Check keywords.
  353.  
  354. miss = 0B
  355. IF (N_Elements(missing) GT 0L) THEN miss = Byte(missing[0])
  356.  
  357. eval_mode = 0
  358. IF (N_Elements(mode) GT 0L) THEN eval_mode = Fix(mode[0])
  359.  
  360. ; *** Set up variables.
  361.  
  362. vx = vol_size1[0]
  363. vy = vol_size1[1]
  364. vz = vol_size1[2]
  365. vxm1 = vx - 1
  366. vym1 = vy - 1
  367. vzm1 = vz - 1
  368.  
  369. ; *** Cell coordinates.
  370. rx = ((vpx2 - vpx1) * Findgen(vx) / Float(vxm1)) + vpx1
  371. ry = ((vpy2 - vpy1) * Findgen(vy) / Float(vym1)) + vpy1
  372. rz = ((vpz2 - vpz1) * Findgen(vz) / Float(vzm1)) + vpz1
  373. xplane = Reform((Temporary(rx) # Replicate(1.0, vy)), (vx * vy))
  374. yplane = Reform((Replicate(1.0, vx) # Temporary(ry)), (vx * vy))
  375. coords = [[Temporary(xplane)], [Temporary(yplane)], $
  376.           [Fltarr((vx * vy), /Nozero)], [Replicate(1.0, (vx * vy))]]
  377.  
  378. ; *** Save current view matrix.
  379. save_t3d = !P.T
  380.  
  381. ; *** Create the volume.
  382. IF (eval_mode GT 0) THEN vol = Bytarr(vx, vy, vz)
  383. IF (eval_mode EQ 0) THEN vol = Intarr(vx, vy, vz)
  384. IF (eval_mode LT 0) THEN vol = Replicate(255B, vx, vy, vz)
  385.  
  386. FOR j=0, (frames-1) DO BEGIN
  387.    ; *** Set up transform.
  388.    T3d, /Reset
  389.    T3d, Rotate=[obj_rot1[0, j], 0.0, 0.0]
  390.    T3d, Rotate=[0.0, obj_rot1[1, j], 0.0]
  391.    T3d, Rotate=[0.0, 0.0, obj_rot1[2, j]]
  392.    T3d, Translate=obj_pos1[*, j]
  393.    IF (focal1[j] GT 0.0) THEN BEGIN
  394.       T3d, Translate=[0.0, 0.0, -dist1[j]]
  395.       T3d, Perspective=$
  396.            ((focal1[j] / obj_pos1[2, j]) * (obj_pos1[2, j] - dist1[j]))
  397.    ENDIF
  398.    T3d, Scale=[img_mag1[0, j], img_mag1[1, j], 1.0]
  399.    T3d, Translate=[img_ref1[0, j], img_ref1[1, j], 0.0]
  400.  
  401.    ; *** Fill in the volume one plane at a time.
  402.    FOR i=0, vzm1 DO BEGIN
  403.       coords[0, 2] = Replicate(rz[i], (vx * vy))
  404.       plane_index = coords # !P.T
  405.       vol_x  = plane_index[*, 0] / plane_index[*, 3]
  406.       vol_y  = plane_index[*, 1] / plane_index[*, 3]
  407.  
  408.       IF (eval_mode LT 0) THEN $
  409.          vol[0, 0, i] = vol[*, *, i] < $
  410.             Interpolate(images[*, *, j], vol_x, vol_y, Missing=miss, $
  411.                         Cubic=Keyword_Set(cubic))
  412.  
  413.       IF (eval_mode GT 0) THEN $
  414.          vol[0, 0, i] = vol[*, *, i] > $
  415.             Interpolate(images[*, *, j], vol_x, vol_y, Missing=miss, $
  416.                         Cubic=Keyword_Set(cubic))
  417.  
  418.       IF (eval_mode EQ 0) THEN $
  419.          vol[0, 0, i] = vol[*, *, i] + $
  420.             Fix(Interpolate(images[*, *, j], vol_x, vol_y, Missing=miss, $
  421.                             Cubic=Keyword_Set(cubic)))
  422.  
  423.    ENDFOR
  424.    Print, 'Completed image ' + String(j+1)
  425. ENDFOR
  426.  
  427. IF (eval_mode EQ 0) THEN vol = Byte(vol / frames)
  428.  
  429. ; *** Restore view matrix.
  430. !P.T = save_t3d
  431.  
  432. RETURN, vol
  433. END
  434.